Chapter 9 Cokriging
9.1 Librerías
library(sp)
library(gstat)
library(sf)
library(rgdal)
library(ggplot2)
library(plotly)##
## Attaching package: 'plotly'
## The following object is masked from 'package:MASS':
##
## select
## The following object is masked from 'package:vcdExtra':
##
## summarise
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
library(Matrix)9.2 Descripción de los datos
Cokriging para las variables \(NO2\), \(O3\), y \(NOX\). La variable de principal riesgo es ozono (\(O3\)), así que se usan las otras dos como covariables espaciales. Día 2020/01/16 A las 17 horas.
datos <- read.csv("2_COK_G_stat/Air_polution_cdmx_2020_01_16_17h.csv")
datos <- datos[c("Estacion",
"X",
"Y",
"NO2",
"O3",
"NOX")]
pander::pander((datos))| Estacion | X | Y | NO2 | O3 | NOX |
|---|---|---|---|---|---|
| AJU | 482901 | 2117907 | NA | 50 | NA |
| AJM | 478188 | 2130946 | 5 | 51 | 8 |
| ATI | 473346 | 2164689 | 20 | 70 | 20 |
| CAM | 482180 | 2152665 | 23 | 83 | 24 |
| CCA | 481502 | 2136931 | 6 | 46 | 8 |
| CUA | 469366 | 2141275 | 9 | 56 | 10 |
| CUT | 479189 | 2180751 | 13 | 75 | 14 |
| FAC | 474444 | 2154232 | 31 | 71 | 35 |
| HGM | 484020 | 2146380 | 25 | 81 | 31 |
| IZT | 487647 | 2143367 | 21 | 61 | 23 |
| LLA | 495842 | 2164872 | 14 | 81 | 15 |
| LPR | 487650 | 2160000 | NA | 64 | NA |
| MER | 487445 | 2147815 | 23 | 79 | 26 |
| MGH | 478716 | 2145543 | NA | 66 | NA |
| MON | 510196 | 2151776 | 9 | 80 | 10 |
| MPA | 498809 | 2123036 | 1 | 50 | NA |
| NEZ | 497038 | 2144394 | 8 | 62 | 9 |
| PED | 478557 | 2136817 | 6 | 56 | 6 |
| SAG | 496819 | 2159801 | 16 | 80 | 16 |
| SFE | 472393 | 2140390 | 9 | 66 | 10 |
| TAH | 498890 | 2128098 | 3 | 47 | 3 |
| TLA | 478535 | 2159383 | 26 | 64 | 30 |
| TLI | 481421 | 2167509 | 18 | 86 | 20 |
| UAX | 489113 | 2134517 | NA | 54 | NA |
| UIZ | 492241 | 2140751 | 7 | 45 | 7 |
| VIF | 489875 | 2173664 | 19 | 80 | 22 |
| XAL | 491355 | 2159031 | 27 | 80 | 27 |
9.3 Matrices de coregionalización.
9.3.1 Matriz definida positiva para el modelo Esférico.
mat1 <- cbind(c(30, 30, 30),
c(30, 50, 30),
c(30, 30, 35))
#matriz definida positiva "cercana"
mat1 <- data.frame(as.matrix(nearPD(mat1)$mat))
names(mat1) <- c("NO2", "O3", "NOX")
row.names(mat1) <- c("NO2", "O3", "NOX")
pander::pander(mat1)| NO2 | O3 | NOX | |
|---|---|---|---|
| NO2 | 30 | 30 | 30 |
| O3 | 30 | 50 | 30 |
| NOX | 30 | 30 | 35 |
9.3.2 Matriz definida positiva para el modelo efecto Hueco.
mat2 <- cbind(c(13.02, 24.5, 18.739),
c(24.58, 46.4, 35.36),
c(18.73, 35.36, 26.95))
mat2 <- data.frame(as.matrix(nearPD(mat2)$mat))
names(mat2) <- c("NO2", "O3", "NOX")
row.names(mat2) <- c("NO2", "O3", "NOX")
pander::pander(mat2)| NO2 | O3 | NOX | |
|---|---|---|---|
| NO2 | 13.02 | 24.54 | 18.73 |
| O3 | 24.54 | 46.4 | 35.36 |
| NOX | 18.73 | 35.36 | 26.96 |
9.4 Definición de objeto en gstat
9.4.1 Semivariogramas univariados
vgmno2 <- vgm(psill = mat1[1, 1],
model = "Sph",
range = 6096,
add.to = vgm(psill = mat2[1, 1],
model = "Hol",
range = 2294))
vgmo3 <- vgm(psill = mat1[2, 2],
model = "Sph",
range = 6096,
add.to = vgm(psill = mat2[2, 2],
model = "Hol",
range = 2294))
vgmnox <- vgm(psill = mat1[3, 3],
model = "Sph",
range = 6096,
add.to = vgm(psill = mat2[3, 3],
model = "Hol",
range = 2294))9.4.2 Semivarogramas cruzados (Bivariados)
vgmno2_o3 <- vgm(psill = mat1[1, 2], model = "Sph",
range = 6096,
add.to = vgm(psill = mat2[1, 2],
model = "Hol",
range = 2294))
vgmno2_nox <- vgm(psill = mat1[1, 3],
model = "Sph",
range = 6096,
add.to = vgm(psill = mat2[1, 3],
model = "Hol",
range = 2294))
vgmno3_nox <- vgm(psill = mat1[2, 3],
model = "Sph",
range = 6096,
add.to = vgm(psill = mat2[2, 3],
model = "Hol",
range = 2294))9.4.3 gstat
remove_na <- function(frame, vari_) {
# Remove na from sp object
datos1 <- frame
bool <- !is.na(datos1@data[vari_])
datos1@data <- datos1@data[bool, ]
datos1@coords <- datos1@coords[bool, ]
return(datos1)
}
coordinates(datos) <- ~ X + Y
g_st <- gstat(NULL,
id = "NO2",
formula = NO2 ~ X + Y,
model = vgmno2,
data = remove_na(datos, "NO2"))
g_st <- gstat(g_st,
id = "O3",
formula = O3 ~ Y,
model = vgmo3,
data = remove_na(datos, "O3"))
g_st <- gstat(g_st,
id = "NOX",
formula = NOX ~ Y,
model = vgmnox,
data = remove_na(datos, "NOX"))
#Cruzados
g_st <- gstat(g_st,
id = c("NO2", "O3"),
model = vgmno2_o3)
g_st <- gstat(g_st,
id = c("NO2", "NOX"),
model = vgmno2_nox)
g_st <- gstat(g_st,
id = c("O3", "NOX"),
model = vgmno3_nox)
pander::pander(do.call(rbind, g_st$model)[, 1:3])| model | psill | range | |
|---|---|---|---|
| NO2.1 | Hol | 13.02 | 2294 |
| NO2.2 | Sph | 30 | 6096 |
| O3.1 | Hol | 46.4 | 2294 |
| O3.2 | Sph | 50 | 6096 |
| NOX.1 | Hol | 26.96 | 2294 |
| NOX.2 | Sph | 35 | 6096 |
| NO2.O3.1 | Hol | 24.54 | 2294 |
| NO2.O3.2 | Sph | 30 | 6096 |
| NO2.NOX.1 | Hol | 18.73 | 2294 |
| NO2.NOX.2 | Sph | 30 | 6096 |
| O3.NOX.1 | Hol | 35.36 | 2294 |
| O3.NOX.2 | Sph | 30 | 6096 |
9.4.4 Estimación del semivariograma
plot(variogram(g_st),
model = g_st$model,
pl = T,
xlab = "Distancias",
ylab = "Semivarianza")
9.4.5 Mapas de predicción de O3 con las covariables espaciales NO2 y NOX
prediction_plot <- function(g_object, variable, map_path) {
map <- readOGR(map_path)
new <- sp::spsample(map, n = 100000, type = "regular")
coordinates(new) ~ x1 + x2
colnames(new@coords) <- c("X", "Y")
predic <- predict(g_object, newdata = new)
prediction <- data.frame(predic)
pred <- paste(variable, ".pred", sep = "")
plot <- ggplot(prediction, aes_string("X", "Y", fill = pred)) +
geom_tile() +
scale_fill_viridis_c() +
theme_void()
return(plot)
}
variance_plot <- function(g_object, variable, map_path) {
map <- readOGR(map_path)
new <- sp::spsample(map, n = 10000, type = "regular")
coordinates(new) ~ x1 + x2
colnames(new@coords) <- c("X", "Y")
predic <- predict(g_object, newdata = new)
prediction <- data.frame(predic)
var <- paste(variable, ".var", sep = "")
plot <- ggplot(prediction, aes_string("X", "Y", fill = var)) +
geom_tile() +
scale_fill_viridis_c(option = "inferno",
direction = -1) +
theme_void()
return(plot)
}
cv_plot <- function(g_object, variable, map_path) {
map <- readOGR(map_path)
new <- sp::spsample(map, n = 10000, type = "regular")
coordinates(new) ~ x1 + x2
colnames(new@coords) <- c("X", "Y")
predic <- predict(g_object, newdata = new)
prediction <- data.frame(predic)
pred <- paste(variable, ".pred", sep = "")
var <- paste(variable, ".var", sep = "")
aux <- abs(sqrt(prediction[var]) / abs(prediction[pred]))
aux[aux > 1] <- 1
prediction["cv"] <- aux
plot <- ggplot(prediction, aes_string("X", "Y", fill = "cv")) +
geom_tile() +
scale_fill_viridis_c(option = "magma",
direction = -1) +
theme_void()
return(plot)
}
pl1 <- prediction_plot(g_st, "O3",
"2_COK_G_stat/SP/mpiosutm.shp")## OGR data source with driver: ESRI Shapefile
## Source: "/home/martha/Documentos/Cursos EE UN/2_COK_G_stat/SP/mpiosutm.shp", layer: "mpiosutm"
## with 54 features
## It has 7 fields
## Linear Model of Coregionalization found. Good.
## [using universal cokriging]
pl2 <- variance_plot(g_st, "O3",
"2_COK_G_stat/SP/mpiosutm.shp")## OGR data source with driver: ESRI Shapefile
## Source: "/home/martha/Documentos/Cursos EE UN/2_COK_G_stat/SP/mpiosutm.shp", layer: "mpiosutm"
## with 54 features
## It has 7 fields
## Linear Model of Coregionalization found. Good.
## [using universal cokriging]
pl3 <- cv_plot(g_st, "O3",
"2_COK_G_stat/SP/mpiosutm.shp")## OGR data source with driver: ESRI Shapefile
## Source: "/home/martha/Documentos/Cursos EE UN/2_COK_G_stat/SP/mpiosutm.shp", layer: "mpiosutm"
## with 54 features
## It has 7 fields
## Linear Model of Coregionalization found. Good.
## [using universal cokriging]
ggplotly(pl1)